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In  this  work  we  estimate  the  state  of  charge  (SOC)  of  Ni-MH  rechargeable  batteries  using  the  Kalman 
filter  based  on  a  simplified  electrochemical  model.  First,  we  derive  the  complete  electrochemical  model 
of  the  battery  which  includes  diffusional  processes  and  kinetic  reactions  in  both  Ni  and  MH  electrodes. 
The  full  model  is  further  reduced  in  a  cascade  of  two  parts,  a  linear  time  invariant  dynamical  sub-model 
followed  by  a  static  nonlinearity.  Both  parts  are  identified  using  the  current  and  potential  measured  at  the 
terminals  of  the  battery  with  a  simple  1-D  minimization  procedure.  The  inverse  of  the  static  nonlinearity 
together  with  a  Kalman  filter  provide  the  SOC  estimation  as  a  linear  estimation  problem.  Experimental 
results  with  commercial  batteries  are  provided  to  illustrate  the  estimation  procedure  and  to  show  the 
performance. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  development  of  low-cost  high  capacity  batteries  is  of  great 
importance  for  large-scale  production  and  application.  Ni-MH  bat¬ 
teries  based  on  a  hydrogen  storage  alloy  have  received  and  still 
receive  much  attention  because  of  their  higher  energy  density, 
superior  charge-discharge  characteristics  and  freedom  from  toxic 
materials.  The  state  of  charge  (SOC)  of  a  rechargeable  battery  is  an 
important  quantity  as  it  is  a  measure  of  the  amount  of  electrical 
energy  available.  In  order  to  guarantee  a  good  performance  and 
lifetime  of  the  battery,  the  SOC  is  usually  kept  within  appropriate 
limits,  for  example  20%  <  SOC  <  90%,  so  the  estimation  of  the  SOC  is 
essential  for  the  battery  to  operate  within  these  safe  limits.  Battery 
SOC  estimation  using  available  potential  and  current  measured  at 
the  terminals  of  the  battery  is  the  goal  of  this  paper. 

There  have  been  several  approaches  to  estimate  the  state  of 
charge  of  a  battery.  In  [1]  a  review  of  impedance  measurements 
to  estimate  SOC  is  presented.  These  approaches  are  based  on  the 
variation  of  the  frequency  response  of  the  battery  at  different  states 
of  charge  and  other  physical  variables.  In  [2]  an  observer  for  linear 
time  variant  systems  was  used  to  estimate  the  open  circuit  volt¬ 
age  (OCV),  which  is  related  to  the  SOC,  from  an  empirical  battery 
model.  Assuming  observability  of  an  extended  model  formulation, 
both  parameters  and  state  of  charge  are  estimated  at  once.  Plett  [3] 
uses  a  simple  dynamical  model  of  one  state,  the  SOC,  which  is  in 
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fact  the  integral  of  the  current.  An  empirical  relationship  between 
SOC  and  measured  current  and  potential  is  used  together  with  the 
extended  Kalman  filter  (EKF).  Further  improvements  were  obtained 
by  adding  new  states  related  to  the  diffusional  process  inside  the 
electrodes.  In  a  different  approach,  Barbarissi  et  al.  [4],  presents 
the  SOC  estimation  based  on  an  electrochemical  model  taking  into 
account  the  kinetics  of  the  reactions  and  the  diffusional  process 
inside  the  electrode.  The  estimator  used  in  this  case  consists  of  a 
simulation  of  the  electrochemical  model  with  the  measured  cur¬ 
rent/potential  as  inputs.  SOC  is  estimated  with  the  internal  variables 
obtained  by  simulation. 

In  this  work  we  estimate  the  SOC  of  Ni-MH  rechargeable  bat¬ 
teries  based  on  an  electrochemical  model  and  using  the  Kalman 
filter.  The  main  contributions  of  this  paper  are  the  following:  First, 
the  derivation  of  a  complete  model  including  both  electrodes,  Ni 
and  MH,  absorption/adsorption  processes  and  double  layer  capac¬ 
ity.  Second,  the  derivation  of  an  electrochemical  reduced  order  of 
the  battery  as  a  Wienner  model  formulation.  It  consists  of  two  parts, 
a  linear  dynamic  in  cascade  with  a  static  nonlinearity.  Third,  we  pro¬ 
pose  an  identification  procedure  by  using  a  simple  1  -D  optimization 
procedure.  Fourth,  the  SOC  estimation  is  performed  in  the  frame¬ 
work  of  linear  systems  in  which  the  estimation  error  bound,  due 
to  the  presence  of  disturbances,  can  be  optimized  using  a  Kalman 
filter.  The  paper  is  organized  as  follows:  In  Section  2,  the  complete 
model  of  both  electrodes  is  presented  including  kinetics  of  reac¬ 
tions,  hydrogen  diffusion  process,  and  charge  balances.  Considering 
the  fact  that  it  will  be  used  for  SOC  estimation,  the  simpler  reduced 
order  model  is  derived  from  the  complete  one.  In  Section  3,  the 
functional  relationship  between  SOC  and  hydrogen  concentration 
is  presented  for  different  electrode  geometries.  In  Section  4,  the 
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Nomenclature 

an 

effective  transfer  area  of  Ni  interface 

dm 

effective  transfer  area  of  MH  interface 

En(t) 

potential  at  Ni-electrolyte  interface 

£m(£) 

potential  at  MH-electrolyte  interface 

Ebat 

potential  at  the  battery  terminals 

Eeq  . 

a,  ref 

equilibrium  potential  of  Ni  interface  at  x„ef(0,  t)  = 

0.5 

Eeq  . 

m,ref 

Equilibrium  potential  of  MH  interface  at  x^f(0,  t)  = 

0.5 

F 

Faraday  constant 

h 

sampling  period 

yo 

Faradaic  current  of  Ni  electrode 

yo 

Faradaic  current  of  MH  electrode 

Icn 

double  layer  capacitive  current  at  Ni  interface 

Icm 

double  layer  capacitive  current  at  MH  interface 

^bat 

current  at  the  battery  terminals 

1 0 

Ire/ 

constant  exchange  current  at  Ni  interface 

7° 

2ref 

constant  exchange  current  at  MH  interface 

Jn{z,  t) 

flux  of  hydrogen  at  position  z  of  Ni  electrode 

Jm{Z,  t) 

flux  of  hydrogen  at  position  z  of  MH  electrode 

Jn (Zj,  t) 

flux  of  hydrogen  at  ith  slice  of  Ni  electrode 

Jm{Zi,  t) 

flux  of  hydrogen  at  ith  slice  of  MH  electrode 

i<i(t ) 

reaction  rate 

L 

order  of  the  Taylor  expansion  series 

t 

time 

Xn (z,  t) 

H  fractional  concentration  at  position  z  of  Ni  elec¬ 
trode 

Xm(.Z,  t ) 

H  fractional  concentration  at  position  z  of  MH  elec¬ 
trode 

Xn(Zi,  t ) 

H  fractional  concentration  at  ith  slice  of  Ni  electrode 

Xm(Zj,  t) 

H  fractional  concentration  at  ith  slice  of  MH  elec¬ 
trode 

*8(0 

short  notation  of  xn(z0,  t) 

Xj(t) 

ith  state  of  simplified  model 

x(t) 

state  vector  of  simplified  model 

x(kh) 

discrete-time  state  vector  of  simplified  model 

x(kh) 

Estimated  state  vector  using  Kalman  filter 

x0  (kh) 

discrete-time  state  x0(t) 

xnef(0,  t ) 

H  reference  fractional  concentration  at  Ni  interface 

*!nef(0,  t) 

H  reference  fractional  concentration  at  MH  interface 

X„(t) 

vector  of  concentrations  in  Ni 

Xm(t) 

vector  of  concentrations  in  MH 

x{kh) 

estimation  error  vector 

z 

spatial  position  along  axis  z 

Zi 

discrete-spatial  position  of  slice  i  along  axis  z 

yo 

incremental  potential  at  Ni-electrolyte  interface 

rjm(t) 

incremental  potential  at  MH-electrolyte  interface 

m 

fractional  surface  concentration  of  MHad 

&ref 

surface  concentration  at  the  equilibrium  reference 
state 

@g 

parameter  vector  of  the  static  nonlinearity 

9gi 

parameter  vector  of  the  software  sensor 

simplified  model  is  identified  using  a  simple  1-D  optimization  pro¬ 
cedure.  The  parameters  of  a  proposed  software  sensor  -  to  estimate 
the  concentration  at  the  surface  of  the  Ni  electrode  to  be  used  in 
SOC  estimation  -  is  also  obtained  in  this  section.  In  Section  5,  a 
linear  Kalman  filtering  strategy  for  the  state  estimation  of  hydro¬ 
gen  concentration  is  presented.  Three  experimental  examples  are 
used  to  illustrate  the  procedure  and  to  show  the  performance  of  the 
proposed  technique. 


2.  Model  formulation 

The  nickel/metal  hydride  battery  consists  of  a  nickel  positive 
electrode  and  a  metal  hydride  negative  electrode  separated  by  a 
porous  inert  separator,  being  the  whole  assembly  immersed  in  a 
highly  concentrated,  usually  30  percent  of  the  weigh,  KOH  aqueous 
electrolyte.  In  Fig.  1  a  scheme  is  shown.  The  main  reactions  during 
discharge  are  the  reduction  of  nickel  oxyhydroxide  in  the  positive 
electrode  and  the  oxidation  of  the  metal  hydride  in  the  negative 
one,  [5].  An  adsorption/absorption  reaction  mechanism  governs  the 
insertion  of  hydrogen  atoms  in  the  metal  [6,7].  The  set  of  reactions 
are  given  as  follows: 

/q  ,disch. 

NiOOH  +  H20  +  e-  =  Ni(OH)2+OH“  (1) 

I<2.disch. 

MHad  +  OH“  H20  +  M  +  e-  (2) 

k®,disch. 

M  +  SHab  MHad  +  S,  (3) 

being  the  reaction  rates  7Q(t)  =  /ci°exp(bo'z£(t))  for  i  =  1, 2  for  elec¬ 
trochemical  steps  in  (1)  and  (2)  with  constants  and  b  =  F/RT , 
is  the  constant  reaction  rate  due  to  the  chemical  reaction  in  (3), 
E(t )  the  potential  at  the  electrode  interface,  being  En(t)  and  Em(  t)  the 
potential  at  Ni  and  MH  interfaces,  respectively,  and  ate[  0,1]  are  the 
symmetry  factors.  MHad  is  an  adsorbed  hydrogen  atom  at  an  active 
site  on  the  alloy  surface,  and  SHab  is  a  hydrogen  atom  absorbed  at 
an  interstitial  site  in  the  alloy,  located  just  beneath  the  surface,  and 
S  is  an  empty  interstitial  site.  The  process  is  reversed  during  charge 
with  rate  constants  K_z(t)  =  /c^exp^cq  -  1  )F(t))  for  i  =  1,2,  and 
/c° 3  for  reactions  (1),  (2),  and  (3),  respectively.  There  are  also  sec¬ 
ondary  reactions  due  to  oxygen  evolution  at  the  nickel  electrode 


Nickel  electrode  MH  electrode 


Fig.  1.  Schematic  diagram  of  the  Ni/MH  cell. 
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and  oxygen  reduction  at  the  metal  hydride  electrode.  However,  in 
the  operative  SOC  range  of  the  battery  (10-90%)  the  secondary  reac¬ 
tions  are  insignificant  compared  with  the  main  reaction.  Since  we 
are  interested  in  modeling  the  dynamics  in  such  range,  the  sec¬ 
ondary  reactions  will  be  neglected  in  this  work.  Thus,  the  model 
considers  the  kinetics  of  the  reactions  taking  place  at  the  electrodes 
described  by  Eqs.  (1  )-(3).  In  order  to  complete  the  dynamics  of  the 
electrochemical  processes  involved,  mass  transport  of  H  atoms  in 
the  metal  hydride  particles  and  Ni  active  material  must  be  taken 
into  account.  Mass  transport  shall  be  modeled  as  a  diffusional  pro¬ 
cess  and  described  in  terms  of  Fields  laws.  We  will  describe  each 
process  in  detail. 

Considering  the  global  process  of  reaction  (1 )  as  the  reaction  of 
a  “free  hydrogen  site”,  NiOOH,  with  H20  to  produce  an  “occupied 
hydrogen  site”,  Ni(OH)2  and  OH-,  the  faradaic  current  related  to 
reaction  (1 )  at  the  Ni  electrode  is  governed  by  the  following  charge 
balance  [8,9]: 


//l(t)  =  F(K1(t)x„(0,t)-IC_1(t)(l  -x„(0,t))),  (4) 

where  xn(0,  t)  =  cn(0,  t)/cn(0,  t)  is  the  fractional  concentration  of 
hydrogen  atoms  in  the  Ni  active  material  at  the  electrochemical 
interface  (Ni  oxide/electrolyte  interface,  named  spatial  position 
z  =  0  at  time  t).  The  hydrogen  concentration,  cn(z,  t),  and  its  maxi¬ 
mum  admissible  value,  cn(z,  t),  define  xn(z,  t)  =  1  and  xn(z,  t)  =  0 
for  maximum  and  zero  concentration  at  position  z  and  time  t, 
respectively.  Let  us  introduce  the  incremental  potential  rjn(t)  = 
En(t)  -  E„qref  as  the  difference  between  the  electrode  potential  and 

the  equilibrium  potential  at  the  reference  state  xj^Zo,  t)  =  0.5. 
Being  the  exchange  current  / ^  defined  by  the  expression 

J°re/  =  Fk^eba'E^ref  =  (5) 

by  replacing  rjn(t)  in  the  expression  of  K\  and  K_\  in  (4)  and  using 
(5)  the  current  at  the  Ni  electrode  can  be  rewritten  as 

Ifi  (t)  =  tfref  (e^“'^n(t)Xn(0,  t)  —  -X„(0 ,!)))•  (6) 

Using  the  same  notation,  the  faradaic  current  balance  at  the  MH 
electrode  is  given  by  reaction  in  (2)  and  can  be  written  as 

J/2(t)  =  -F  (K2(tWt)  -  K_2(t)(  1  -  0(t)))  ,  (7) 

being  0(t)  the  fractional  surface  concentration  of  MHad  species.  As  in 
the  Ni-electrode,  we  define  rjm(t)  =  Em(t)  -  E6^^  as  the  difference 
between  the  electrode  potential  and  the  equilibrium  potential  at 
the  reference  state  x^f(0,  t)  =  0.5  (the  fractional  concentration  of 
hydrogen  atoms,  SHab,  at  the  metal  hydride  particles-electrolyte 
interface).  Accordingly,  the  exchange  current  is  given  by 

IC2 lef  =  Fk°eb0l2E™.referef  =  FkVb(“2_1)^re/(l  -  Ortf),  (8) 

where  0ref  is  the  surface  concentration  at  the  equilibrium  reference 
state.  Taking  into  account  the  fast  dynamical  response  of  the  hydro¬ 
gen  absorption  reaction  (HAR  step  given  by  reaction  3),  equilibrium 
may  be  assumed  and  the  following  relationship  is  fulfilled  [10]: 

k°(l  -0(t)]xm(O,t)  =  fc°30(t)(l  -xm(0 ,t)).  (9) 


By  replacing  0(t)  from  (9)  in  (7)  and  taking  into  account  that 
0ref  =  1  /{K  +  1 )  with  I<  =  k°3//<3  -  obtained  with  x^f(0,  t)  =  0.5  in 
(9)-  and  using  (8),  the  expression  of  the  faradaic  current  at  the  MH 
electrode  can  be  written  as 


I  _  _ 2, re/ _ 

h  xm(z0,t)  +  K(\ -xm(z0,t)) 


(xm(0,  t)eb0l^m 


(10) 


Fig.  2.  Spatial  discretization  of  the  active  film. 


Considering  the  current  during  discharge  to  be  positive,  the 
complete  current  and  potential  balance  is  given  by  the  series  of 
both  electrodes  fulfilling 


hat  ~  ~hi  ~  hn  —  If2  +  I  cm  (H  ) 

Ebat  —  En  —  Em  —  hatE-i •  (12) 


All  the  variables  are  function  of  time.  Ibat  and  Ebat  stand  for  the 
current  and  potential  at  the  battery  terminals;  Rt  is  the  resistance  of 
the  electrolyte;  and  /cn,  Icm  are  capacitive  currents  due  to  the  double 
layer  at  the  Ni  and  MH  electrode  interface  respectively,  having  the 
following  expression: 


—  dEn  _  clEm 

lcn  ~  ~dFLdl,nt  lcm  ~  ~dFLdl,mi 


(13) 


where  Q/  n  and  Q/  m  are  the  double  layer  capacities  of  the  electrode 
surfaces.  The  dynamics  of  hydrogen  concentration  in  the  electrodes 
are  given  by  diffusional  processes  described  by  Fields  laws.  The  Ni 
electrode,  represented  schematically  in  Fig.  1 ,  shall  be  described  as  a 
porous  metallic  structure,  where  the  active  material  is  deposited  on 
the  cylindrical  pores  of  mean  dimensions  [5].  The  active  material  on 
the  pore  walls  forms  a  cylindrical  film  of  average  thickness,  d  =  Re  - 
Rif  hydrogen  diffusion  taking  place  mainly  in  the  radial  direction. 
In  a  similar  way,  but  considering  spherical  geometry  particles,  the 
hydrogen  diffusion  in  metal  hydride  electrode  is  governed  by  Fields 
laws.  The  analytical  solution  of  Fields  laws  is  complex  [11  ].  Then,  in 
order  to  have  a  working  expression  of  the  diffusional  processes,  they 
can  be  approximated  into  a  set  of  ordinary  differential  equations 
by  using  a  spatial  discretization.  The  spatial  discretization  is  a  very 
well-known  method  to  approximate  partial  differential  equations 
in  ordinary  differential  equations.  For  details  see  [12,13],  and  also 
applied  to  Ni-MH  batteries  in  [4].  Each  diffusional  process  can  be 
discretized  along  the  space  variable  z  by  considering  N  slices  of  the 
metal  with  thickness  A z,  as  illustrated  in  Fig.  2.  By  denoting  spatial 
position  of  ith  slice  as  Z/,  the  dynamics  of  diffusional  processes  for 
both  electrodes  can  be  written  as  follows  (see  the  Appendix  A  for  a 
detailed  derivation): 


^l=Anxn(0  +  Bn//1(0 
*8(0  —  Enxn(t) 


where  xg(t)  is  a  short  notation  of  xn(z0,  t),  and 


-d0  d0  0 
1  -(1  +di)  d\ 


0  0  0  ■■■  1  — (1+djv-i)  djv_i 

0  0  0-0  1  -1 


(14) 

(15) 


C„  =  [1 0 ■  •  •  0],  B„  =  [F, o,  ■••,0]T;  x„(t)  =  [x„(z0,  t), x„(z, ,  t), . . . ,x„(zN,  t)]r, 

(16) 


(1  -xm(0,  t))eb(“2_1),?m). 
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Fig.  3.  Ni-MH  battery  model. 


with  a  =  Dn/ Az2,  Dn  is  the  diffusion  coefficient,  Az  =  being 
Re  and  Rt  the  exterior  e  interior  radium  in  the  cylindrical  geome¬ 
try,  j3  =  \/{FancnAz)  with  an  the  effective  transfer  area,  and  d,-  = 
Zj+1/Zj  =  (2i  +  3)/(2i  +  1)  for  1  <i<N. 

A  similar  approach  can  be  followed  for  the  metal  electrode  (see 
the  Appendix  A)  by  replacing  Az  =  with  Rp  the  radium  of  the 
sphere,  and  d{  =  z?+1/z?  =  (2 i  +  3)2/(2 i  +  1  )2  for  1  <  i  <  N. 

Finally,  the  complete  battery  model  is  then  given  by  Eqs.  (6), 
(10)-(16).  The  equations  of  the  diffusional  process  for  the  metal 
electrode,  which  are  similar  to  (14)-(16),  are  also  included  in  the 
complete  battery  model  and  depicted  in  Fig.  3. 

2.1  Reduced  complexity  model 

There  are  different  reasons  that  justify  the  use  of  a  simplified 
model  to  be  used  for  SOC  estimation.  We  present  two  important 
arguments:  (i)  In  commercial  Ni-MH  batteries  usually  the  nega¬ 
tive  electrode  (MH  electrode)  has  a  higher  effective  capacity  than 
the  positive  one  (nickel  oxyhydroxide/  nickel  hydroxide  electrode), 
being  quite  reasonable  to  assume  that  variations  of  the  fractional 
hydrogen  concentration  in  the  metal  hydride  electrode  may  be 
neglected  with  respect  to  variations  in  the  Ni  electrode.  Thus,  the 
Ni  electrode  is  sufficient  to  describe  the  dynamics  of  the  energy 
storage  in  the  battery.  Consequently,  the  changes  in  Ebat  are  mainly 
due  to  changes  in  En,  as  reported  in  [5].  Accordingly,  next,  poten¬ 
tial  changes  at  the  MH  electrode  shall  be  neglected,  considering 
Em  =Emref  (ii)  The  maximum  current  is  bounded  for  construc¬ 
tive  reasons.  Then,  charge  and  discharge  are  mainly  due  to  the 
low  frequency  components  of  the  current,  Ibat.  Since  the  current 
passing  throughout  the  double  layer  capacity  contains  mainly  high 
frequency  components,  if  Ibat  is  conveniently  “low  pass  filtered”, 
as  described  in  Section  5.1,  there  are  no  significant  differences  in 
SOC  estimation  by  neglecting  the  capacitive  current  Icn  and  con¬ 


sidering  Ibat(t )  ~  -f/i(0  in  (11).  Then,  in  order  to  get  a  simplified 
model,  we  can  use  Ibat{t)  =  -f/i(t)  in  Eq.  (6)  as  input  and  potential 
Ebat  =  Eenqref  -  Ee^ref  +  r/„(f)  -  lha[R,  as  output.  However,  to  use  this 
approach  it  is  necessary  for  pn{t)  in  (6)  to  be  a  function  of  Xo  and 
Ibat.  In  other  words,  there  should  be  a  unique  value  of  rjn{t)  for  each 
couple  of  values  (xg,  Ibat).  In  order  to  show  that  this  requirement  is 
fulfilled,  let  us  start  by  noting  that  in  (6)  the  concentration  xg  is  a 
function  of  the  couple  ( r\n ,  lbat ),  as  follows: 

„  +  fO  p5[oi\  — l)?7n 

*0  -  eb(aA-\)t]n  +  eba^n  ~  f ^n’  bat^'  ^ 7 ) 

A  continuous  function /is  one  to  one  (and  hence  invertible)  if  and 
only  if  it  is  either  strictly  increasing  or  decreasing  with  not  local 
maxima  or  minima.  In  our  case  the  derivative  with  respect  to  r\n 
after  some  algebraical  manipulation  can  be  written  as 

dxg  b[a i  -  1  )eb“i  **«  ( 1  -  xg )  -  i  be01 1 h ^ 

drjn  e^l  rjn  -\-  eb(a\  ~^)rjn 

Taking  into  account  that  xge(0, 1),  and  aqe(0, 1),  the  deriva¬ 
tive  above  is  always  negative.  Thus,  the  inverse  function  pn  = 
/-1  (x0,  Ibat)  exists,  hence  there  is  only  one  value  of  rjn  for  each 
couple  (xg,/bflt). 

Taking  into  account  these  considerations,  the  complete  model 
depicted  in  Fig.  3  can  be  simplified  as  shown  in  Fig.  4.  We  will  show 
it  later  that  experimental  results  are  accurately  represented  by  the 
proposed  model. 

Remark  1.  The  simplified  model  is  a  cascade  of  a  linear  time 
invariant  dynamical  system  and  a  static  nonlinearity.  Such  kind  of 
model  is  called  Wiener  model.  A  Wiener  model  consists  of  a  linear 
dynamic  system  followed  by  a  static  nonlinearity.  The  input  and 
output  are  measured,  but  not  the  intermediate  signal.  Then,  the 
identification  of  both  parts  as  black  boxes  are  difficult  to  obtain  [14]. 
Using  the  electrochemical  model  we  propose  a  successful  approach 
in  the  next  section. 

Before  to  proceeding  with  the  identification  of  the  model  param¬ 
eters,  we  show  in  next  section  that  the  gain  j3  in  (16)  can  be  known 
directly  from  the  structure  of  the  diffusion  model.  Also  we  show 
how  it  is  related  to  a  working  expression  to  compute  the  SOC. 

3.  The  state  of  charge 

The  SOC  in  the  battery  is  determined  by  integrating  the  current 
as  follows: 

SOC(t)=  — (0o  + 

where  Qo,  and  Qmax  are  the  initial  and  maximum  electrode  charge. 
The  quantity  Qmax  is  called  capacity  of  the  battery  and  it  is  defined  as 
the  product  of  the  constant  current  and  the  time  required  to  reach 
the  maximum  charge,  starting  with  battery  completely  discharged, 
Qo  =  0.  Thus,  SOC  is  the  unity  for  battery  completely  charged  and 
zero  for  completely  discharged.  We  derive  the  expression  for  SOC 
in  the  nickel  electrode  by  integrating,  with  respect  to  time,  both 


haMdr^j 


(19) 


Fig.  4.  Simplified  battery  model. 
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sides  of  equation  (14)  in  which  we  consider  Ibat(t)  =  Since 

in  the  simplified  model  only  the  Ni  electrode  will  be  considered,  the 
subindex  of  the  fractional  concentration  will  be  dropped  for  clarity 
as  follows:  xn(Zj,  t)  will  be  denoted  as  Xj(t),  and  vector  xn(t)  as  x(t). 
After  integrating  and  using  (19),  we  have 


x(t)-x(0) 


Ax(r)dr  -  BQmax(SOC(t)  -  SOC(O)). 


(20) 


Using  the  elements  of  matrix  A  defined  in  (16),  the  previous  is 
equivalent  to  the  following  set  of  equations: 


x0(t)  =  ad0  f  (xi(r)-xo(r))dr-y0Qmax(SOC(t)-SOC(O))  +  xo(O) 
Jo 

xi {t)  =  a  [  (x0(r)-xi(r))  +  di(x2(r)-xi(r))dr  +  xi(0) 


)  =  a  [  ( 
Jo 

i  [  ( 
Jo 


(21) 


XN-1  (t)=a  /  (xN_2(r)-xN_i(r))+div_i(xiv(r)-xiv_i(r))dr+xiv_i(0) 

Jo 

Xjv(0  =  OL  /  (xjv-i(t) - xjv(r))dr. 


By  replacing  the  integral  term  of  xN(t)  from  the  last  equality  in 
the  expression  of  Xjv-i  (t)  we  have 


XN-l(0  =  M 


l 


t 

(xN_2(r)  -  xN_!  (r))dr  -  dN_-,xN(t )  +  xN_!  (0). 


By  following  the  procedure  from  bottom  to  top  on  the  set  of  equa¬ 
tions  (21),  the  following  equality  is  obtained: 


SOC(t)  =  cnt 


xp  (t)  +  dpXi(f)  +  dpdix2(f)  +  •  •  •  +  d0-  •  •djV_iXN(t) 
fiQmax 


(22) 


where  the  constant  cnt  is  given  by 


*o(Q)  +  doXi(O)  +  dpdiX2(0)  +  •  •  •  +  dp-  •  djy_iXjv(0)  +  $oc(0) 
fiQmax 


(23) 


Considering  the  stationary  case,  zero  current  in  the  limit  when  time 
goes  to  infinity,  the  concentrations  x^oo)  for  0  <  i  <  N  are  equal.  For 
battery  completely  charged  all  the  concentrations  are  equal  to  zero 
which  corresponds  to  SOC  equal  to  the  unity  and  vice  versa.  Then, 
from  (22)  a  set  of  two  equations  with  two  unknowns,  cnt  and  /?,  can 
be  formed.  The  solution  gives  cnt  =  1  and 


which  -  using  the  fact  that  SOC  =  1  -  x0(oo)  -  leads  to  the  fulfill¬ 
ment  of  the  following  SOC/OCV  relationship: 


B  —  Ebat 


_  peq  ,  Feq 
m,ref  ^  n,ref 


1 

b 


1  Xq ( OO )  \ 
X0(oo)  ) 


1 

b 


SOC(oo) 

1  -  SOC(oo) 


)  =f-\x o,0) 


(27) 


In  the  next  section  we  identify  the  parameters  of  the  model  to 
be  used  for  SOC  estimation. 

4.  Model  parameter  identification 

The  model  parameters  that  remain  to  be  identified  are  a  in 
(16)  and  the  parameters  corresponding  to  the  static  nonlinearity. 
The  static  nonlinearity  is  a  soft  function  that  can  be  well  approxi¬ 
mated  by  using  different  approaches.  We  choose  to  approximate  the 
function  with  a  Taylor  series  expansion  of  multivariable  functions, 
which  can  be  written  as  a  linear  regression  problem  as  follows: 

Ebat  —  ~^rn,ref  +  ^ n,ref  ~  hat^i  +  /  hat)  —  hat)  — 

(28) 

where  0g  is  a  column  vector  of  parameters  and  cpg ,  a  row  vector  of 
signals  as  follows: 


%  =  [1.  I  bat-  Xo,  Ibat*  o,  liar  *o  kat,  *0’  higher  orders]; 


where  L  depends  on  the  series  order  approximation  of  the  static 
nonlinearity.  In  expression  (29)  the  ith  element  of  0g  is  the  deriva¬ 
tive  of  ith  order  of  the  unknown  function  g(x0,  hat)  evaluated  at 
a  given  coordinates.  We  choose  to  expand  around  the  stationary 
coordinates  (x0,  hat )  =  0).  Using  this  representation,  we  propose 

the  following  procedure  to  estimate  a  and  0g :  (i)  Discharge  com¬ 
pletely  and  recharge  with  a  known  small  amplitude  current  during 
a  given  time.  After  a  long  resting  time,  the  known  stationary  values 
x(0)  are  reached,  (ii)  Excite  the  battery  with  a  time  varying  cur¬ 
rent  and  obtain  a  record  of  M  sampled  of  both  Ibat(kh),  and  Ebat(kh). 
(iii)  Choose  an  arbitrary  value  of  the  scalar  a  within  a  given  griding 
interval,  (iv)  Using  this  value  of  a  in  Eq.  (14),  obtain  the  zero-hold 
discrete-time  system  as  follows: 

x((k  +  \)h)  =  Adx(kh)  -  B dIbatm  (30) 

x0  {kh)  =  Cdx{kh),  (31) 


1  +  do  +  d0di  +  •  •  •  +  do-  •  djv-i 
Qmax 


where  h  is  the  sampling  period,  and  (Ad,  Bd,Cd )  is  the  zero-hold 
discrete-time  system  obtained  from  (An,  Bn,Cn )  as  follows  [15]: 


Finally,  the  expression  for  SOC(t)  is  given  by 


SOC(t)  =  1  -  *0(0  +  do*i(t)  +  d0diX2(t)  +  •  •  •  +  dp-  •  -dN_iXN(t) 

fiQmax 


(25) 


with  given  by  (24). 

Remark  2.  Considering  different  steady  states,  { hat  =  0)  at  the 
equilibrium  (t  oo),  the  charge  balance  (6)  gives  the  relationship 
between  OCV  and  SOC,  often  used  to  compute  the  SOC,  as  follows: 


rh 

Ad  =  eAnh;  Bd=  /  e*s Bnds ;  Cd  =  C„;  (32) 

Jo 

and  x(kh)  is  the  sampled  vector  of  concentrations  at  times  t  =  kh. 
(v)  Using  the  discrete-time  model,  the  M  samples  of  Ibat(kh)>  and 
Ebat(kh),  and  the  known  initial  concentrations,  obtain  the  estimated 
sequence  x0(kh )  by  simulation  using  (30)  and  (31 ).  (vi)  Using  the  M 
values  of  hatU<h)>  Bbat{kh),  and  x0{kh)  obtained  in  the  previous  step, 
the  value  of  parameters  0g  are  fitted  by  using  the  least  squared 
approximation  as  follows: 


x0(oo)eb“1,?n  =  (1  -x0(oo))e'’(“i-1>'" 


(26)  6g  =  (3>t<I>)  1<try, 


(33) 
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where 

$  =  [<PgU)T>  •••>  <Pg(k)T ,  (pg{M)T,]T  (34) 

<Pg(k)=[  1,  Itat(kh),  x0(kh),  lbac(kh)x0(kh), 

I bat (^)»  *o(kh)Ibat(kh),  Xg(fch),  higher  order  elements]  (35) 

y  =  [£tat(/i),  . . . ;  £to(Wi),  . . . ;  Eta(Mfl)]7'.  (36) 

(vii)  Compute  the  sum  of  the  least  squared  error  of  the  estimation 
as 

M 

J2<$  =  (Y-<Mg)T(Y-<Mg).  (37) 

1=1 

Save  the  result,  choose  a  new  value  of  a  and  return  to  step  (iii).  The 
procedure  continues  for  all  the  griding  values  of  a  in  the  consid¬ 
ered  interval.  Finally,  we  adopt  the  coupled  a  and  0g  which  gives 
the  lowest  least  squared  error  subjected  to  the  constraints  that  the 
simulated  concentration  x0(kh)  lies  in  the  interval  [0, 1  ]. 

Remark  3. 

•  It  is  important  to  note  that  the  computational  complexity  is  1- 
D  which  is  solved  by  performing  steps  (iii)-(vii)  repeatedly  for 
different  values  of  a  within  the  grid  interval  of  one  dimension.  It 
must  also  be  noted  that  the  reduced  complexity  is  due  to  the  fact 
that  we  have  used  an  analytical  model  description  of  the  battery. 

•  The  surface  concentration  sampled  sequence,  x0(kh),  is  obtained 
by  simulation  and  it  is  unique  and  different  for  different  values  of 
a.  The  interval  of  possible  values  of  a  can  be  known  from  the  ana¬ 
lytical  expression  a  =  Dn(N  +  1  )2/(fte  -  ft;)2  where  typical  values 
of  Dn  are  in  the  interval  (10-9  to  10-1°),  N  is  given  by  the  number 
of  slices  of  the  spatial  discretization,  and  (fte  -  ft,)  depends  on  the 
battery. 

•  The  identifiability  is  governed  by  the  minimum  least  squared 
problem  (33)  which  states  that  the  estimated  vector  0g  converges 
to  the  true  value  if  the  input  Ibat  is  persistently  excited  (PE),  which 
means  the  current  variations  need  to  be  rich  enough  to  satisfac¬ 
tory  identify  the  parameters  0gl  [16,17]. 

•  Since  the  function  x0  =f(rj,Ibat )  has  inverse  pn  =  /_1(x0, Ibat)* 
then,  from  (28)  the  function  Ebat  =g(xo,Ibat )  has  also  inverse 
given  by  x0  =  g_1(£bat>  hat)-  Then,  the  value  of  parameters  that 
minimizes  the  Euclidean  norm  of  the  cost  in  (37)  also  minimizes 
the  Euclidean  norm  of  the  differences  between  the  real  and  esti¬ 
mated  concentration  x0{kh). 

•  Since  we  are  interested  in  low  frequencies,  a  previous  low  pass 
(LPF)  filtering  of  sequences  Ibat{kh)  and  Ebat(kh)  must  be  done. 
Some  guidelines  about  how  to  choose  the  cut  frequency  of  the 
LPF  will  be  given  in  the  next  section. 

It  is  useful,  for  further  use  in  SOC  estimation,  to  have  a  practical 
way  to  estimate  the  value  ofxo {kh)  from  measurements  Ebat(kh)  and 
Ibatikh).  To  this  end,  we  need  to  identify  the  inverse  function  x0  = 
g -1  (Ebat ,  Ibat )  which  works  in  fact  as  a  software  sensor  of  concentra¬ 
tion  x0{kh).  It  can  be  done  by  using  the  estimated  sequence  x0{kh), 
obtained  in  the  optimization  procedure  described  before,  and  find¬ 
ing  the  coefficients  of  the  linear  regression  model  of  the  inverse 
function  x0{kh)  =  g~\Ebat(kh),  Ibat(kh ))  «=  <Pgi(kh)0gi  =  x™(kh).  The 
estimated  parameters  6gi ,  in  the  sense  of  minimum  least  squared, 
can  be  obtained  in  a  similar  procedure  as  the  one  explained  above 
but  replacing  O  and  Y  in  (33),  by 

<h  =  [^(l)T,  cpgi(kf,  (pgj(M)T ,  f  (38) 


Fig.  5.  SOC  estimation. 


<Pgi(k)  =  [  1,  Ibat(kh),  Ebat(kh\  Ibat(kh)Ebat(kh),  Ibat(kh), 
Ebat{kh)Ibat(kh),  Ebat(kh),  higher  order  elements]  (39) 

y  =  [x0(h),  . . . ,  x0(kh),  . . . ,  xo(Mh)  ]T.  (40) 

Note  that  the  estimated  value  xfi(kh)  =  (pgi{kh)Ogi,  obtained  from 
the  measured  current  4at(/<ft)  and  potential  Ebat(kh ),  is  an  approxi¬ 
mation  of  the  true  value  x0{kh).  It  will  be  used,  in  next  section,  as  a 
software  sensor  of  fractional  concentration  at  the  interface. 

5.  SOC  estimation  and  results 

The  basic  idea  to  estimate  the  SOC  is  to  use  the  expression  (25) 
together  with  (24)  employing  estimated  values  of  the  vector  con¬ 
centrations  x(kh)  that  will  be  named  x(kh).  To  this  end,  let  us  have 
first  the  estimated  value  of  x™{kh)  as  it  was  explained  in  the  last 
section.  It  is  interpreted  as  the  true  value,  Xo {kh),  corrupted  by 
measurement  zero  mean  random  noise,  errors  due  to  the  imper¬ 
fect  interpolation,  and  possible  unmodeled  dynamics.  Using  this 
measure  the  concentration  vector  x{kh)  is  estimated  by  means  of 
the  Kalman  filter.  The  equation  of  the  Kalman  filter  are  given  by 

x((k  + 1  )h)  =  AjX(Wi)  +  BdIbat(kh )  +  Kk(x™(kh)  -  Cdx{kh ))  (41 ) 

Kk=AdPkCd(R2  +  CdPkCTdf 1  (42) 

Pk+ 1  =  AdPkATd  +  R,  -  AdPkCT(R2  +  CdPkCTd  r1  CdPkATd  (43) 

where  P^,  fti ,  and  ft2  are  the  positive  definite  covariance  matrices  of 
errors  in  the  estimations,  noise  disturbances  at  the  concentrations, 
and  measurement  noise,  respectively.  Notice  that  the  global  conver¬ 
gence  is  guaranteed  since  the  diffusional  process  is  linear.  Moreover, 
in  the  ideal  case  where  the  disturbances  at  the  states  and  measure¬ 
ment  noise  can  be  modeled  as  gaussian  white  noise  with  zero  mean, 
the  error  covariance  matrix  between  real  and  estimated  vector  con¬ 
centrations  is  given  by  £(x{kh)  -  x(kh))(x(kh)  -  x{kh))T  =  Pk  in  (43) 
where  S  means  expected  value.  No  more  details  about  Kalman  filter 
will  be  given  since  there  is  a  wide  reference  to  this  respect.  See  for 
example  [15,3].  The  complete  scheme  is  depicted  in  Fig.  5. 

The  measure  xg1  (kh),  provided  by  the  software  sensor  and  based 
on  data  interpolation,  is  an  approximation  of  the  true  state  x0(kh). 
This  error  affects  the  estimation  of  the  vector  x(kh)  in  (41 ).  In  order 
to  analyze  how  the  estimation  is  affected  by  this  error,  let  us  assume 
that  the  measure  x™(kh )  consists  of  the  sum  of  the  true  value,  given 
by  x0{kh)  =  Cdx(kh ),  plus  a  difference  given  by  Ax(/di),  which  is 
in  fact  the  error  introduced  by  the  interpolation  procedure.  Then, 
the  error  between  the  real  and  estimated  concentrations  can  be 
obtained  by  substracting  (41 )  from  (30)  as  follows: 

x{kh  +  h)  =  x(kh  +  h)-  x{kh  +  h) 

=  (Ad-I<kCd)x(kh)  +  KkAx(kh) 


(44) 

(45) 
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Fig.  6.  True  and  low  pass  filtered  current  and  potential  for  experiment  1. 


By  denoting  F  =Ad-  I<kCd,  the  solution  of  the  above  difference 
equation  is  given  by 

k- 1 

x(kh)  =  Fkx{  0)  +  yV-'-^Ax  (46) 

j= o 

It  consists  of  two  terms,  the  first  is  the  transient  due  to  the  error  at 
the  initial  conditions  x(0),  and  the  second  one  is  due  to  the  uncer¬ 
tainties  in  the  software  sensor  Ax.  The  transient  is  governed  by 
the  eigenvalues  of  matrix  F  which  are  tuned  by  the  gain  Kalman 
matrix  I<k  in  a  way  that  short  transient  need  large  values  of  I<k.  The 
second  term  is  a  discrete-time  convolution  between  the  transient 
and  the  product  of  the  gain  matrix  I<k  by  the  interpolation  error 
A  x(/<fr).  Then,  large  values  of  the  gain  matrix  I<k  produces  short  tran¬ 
sient  but  increases  the  effect  of  measurement  errors.  Conversely,  a 
low  filter  gain  produces  a  long  transition  but  mitigates  the  effect  of 
the  uncertainties.  Using  the  covariance  matrix  of  the  measurement 
error  R2  =  £(  AxAj),  the  gain  of  the  Kalman  filter  solve  optimally 


the  tradeoff  between  both  terms  in  (46)  and  it  depends  inversely 
on  R2  according  with  Eq.  (42).  In  the  next  section  we  will  use  this 
analysis  to  interpret  the  experimental  results. 

5.1.  Experimental  results 

The  examples  presented  in  this  paper  are  based  on  measure¬ 
ments  made  in  rechargeable  batteries  Duracell  AAA/HR03/DX2400, 
NiMH/1,2/800  mAh.  The  used  experimental  setup  consists  in  a 
power  unit  driven  by  a  personal  computer  throughout  an  I/O  card 
with  A/D  and  D/A  conversions  in  a  MATLAB  environment.  Three  dif¬ 
ferent  experiments  exciting  the  battery  with  a  variant  current  load 
(charge  and  discharge)  using  a  sampling  time  of  h=  8  s  are  pre¬ 
sented.  The  first  two  experiments  are  used  to  identify  the  model 
parameters  (a,  0g)  together  with  the  software  sensor  parameters 
Ogi.  Both  experiments  are  also  used  to  estimate  SOC.  In  the  third 
experiment,  the  SOC  is  estimated  using  the  values  of  the  model 
parameters  identified  with  the  first  two  experiments. 


Fig.  7.  SOC  estimation  for  experiment  1.  Upper,  with  R2  =  10.  Lower,  with  R2  =  0.1. 
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Fig.  8.  True  and  low  pass  filtered  current  and  potential  for  experiment  2. 


Before  each  experiment  the  battery  was  discharged  completely 
and  recharged  with  a  small  amplitude  known  current.  After  a  long 
resting  time,  the  stationary  values  x(0)  were  perfectly  known.  Then, 
by  integrating  the  current  we  also  knew  the  real  SOC.  For  each 
experiment  current  (Ibat)  and  potential  ( Ebat )  at  the  terminal  of  the 
battery  were  recorded.  In  order  to  estimate  the  SOC,  first  both  cur¬ 
rent  and  potential  samples  were  low  pass  pre-filtered  in  order  to 
eliminate  the  high  frequency  components.  The  filter  used  is  a  zero 
phase  shift  filter  with  the  following  structure: 

100 

sfm  =- T-  ^  s{kh),  (47) 

i=-100 

where  s(kh)  is  the  sampled  signal  and  / (kh)  is  the  filtered  signal. 
This  filtering  task  can  be  done  on  line,  at  the  same  time  that  the  sam¬ 
ples  are  obtained  and  considering,  in  our  case,  the  necessary  delay 
of  100  samples  due  to  the  noncausality  of  the  filter.  The  criterion 
for  choosing  the  cutoff  frequency  of  the  LPF  depends  on  the  battery. 
Using  the  recorded  files  of  current  and  potential,  in  the  identifica¬ 


tion  stage,  we  started  with  a  low  cutoff  frequency  and  iteratively  it 
was  increased  until  the  error  between  the  estimated  and  real  SOC 
reached  some  admissible  bound. 

After  filtering  the  signals  were  re-sampled  at  period  T  =  5.33 
min.  The  current  and  potentials  together  with  their  filtered  version 
for  the  three  experiments  are  shown  in  Figs.  6,  8,  and  10.  Using 
a  partition  number  N  =  30,  the  first  and  second  experiments  were 
used  to  identify  parameter  a.  By  taking  a  value  of  Re  -  Ri  =  10-3  cm, 
the  search  interval  of  the  1-D  optimization  procedure  was  [0.1, 1] 
giving  a  =  0.397  and  its  corresponding  0g  and  0gi.  The  value  of  for 
the  chosen  model  structure  is  2.0855  x  10-4  coulomb. 

With  the  identified  values,  SOC  estimations  were  performed 
using  the  first  and  second  experiments.  In  all  the  experiments 
the  initial  concentration  values  of  the  Kalman  filter  are  arbitrar¬ 
ily  fixed  in  0.1,  0.5,  and  0.9.  There  were  used  two  different  values 
of  matrix  covariance  of  the  measurement  error,  R2  =  10,  and  R2  = 
0.1,  for  computation  of  Kalman  gain  matrix  I<k.  Large  value  of 
R2  gives  small  filter  gain  with  relatively  long  transient  and  rel¬ 
atively  good  match  to  the  true  SOC.  Opposite,  small  value  of  R2 
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gives  high  gain  value  with  short  transient  but  with  a  poor  fit  of 
the  real  SOC,  due  to  the  amplification  of  the  measurement  error. 
In  Figs.  7  and  9  the  estimated  SOC,  using  both  gain  values,  are 
depicted  for  first  and  second  experiments.  The  optimal  gain  value 
can  be  obtained  by  knowing  the  value  of  R2.  It  can  be  obtained 
from  the  interpolation  procedure  in  the  identification  stage.  It  must 
be  noted  that  as  the  interpolation  error  is  smaller,  the  better  the 
estimates. 

The  study  of  different  interpolation  methods  than  the  linear 
regression  used  in  this  work,  like  splines  [19],  kernel  [18],  or  Krig- 
ing  [20]  between  others,  to  obtain  software  sensors  with  minimum 
error  will  be  grounds  for  future  research. 

With  the  model  parameters  and  smoothing  filter  properly 
selected  in  the  identification  stage,  the  SOC  in  third  experiment 
is  performed  and  depicted  in  Fig  11.  Also  two  different  gain  matrix, 
R2  =  10,  and  R2  =  0.1,  was  employed  in  order  to  show  the  effect  on 
the  estimates.  It  must  be  noted  that  a  good  fitting  of  software  sensor 
parameter  is  a  crucial  point  to  improve  the  results. 


6.  Conclusions 

A  complete  electrochemical  model  of  two  electrodes  of  Ni-MH 
rechargeable  batteries  for  SOC  estimation  was  presented.  Taking 
into  account  constructive  reasons,  the  two-electrode  model  has 
been  reduced  to  a  simpler  one-electrode  model  which  belongs  to 
the  family  of  the  Wiener  models.  The  simplified  model  takes  into 
account  only  the  slow  current  variations,  which  are  in  fact  the  com¬ 
ponents  that  change  the  SOC.  Using  the  kinetic  equations,  the  two 
parts  of  the  Wiener  model,  the  dynamic  linear  and  the  static  nonlin¬ 
earity  were  identified  by  using  1  -D  optimization  method  to  tune  the 
parameters.  The  interval  search  is  bounded  by  previous  knowledge 
of  typical  parameter  values. 

The  proposed  method  allows  the  estimation  of  the  state  of 
charge  in  Ni-MH  rechargeable  batteries  using  a  simplified  electro¬ 
chemical  based  model  by  means  of  a  software  sensor  and  a  Kalman 
filter  to  estimate  the  concentration  inside  the  Ni  electrode.  Using 
this  setup,  minimum  covariances  of  the  estimation  errors  between 
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real  and  estimated  concentrations  due  to  external  disturbances  are 
obtained.  Also  the  bound  of  the  estimation  errors  are  known.  The 
expression  of  SOC  as  a  function  of  the  estimated  concentrations 
for  different  electrode  particles  geometry  was  given.  The  method 
was  used  in  experiments  with  commercial  batteries  reaching  the 
expected  performance. 
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Appendix  A. 


The  diffusion  process  in  the  Ni  electrode  is  described  by  the 
cylindrical  form  of  Fields  laws  considering  only  the  radial  direction 
(coordinate  z)  as  follows: 


dx„(z,  t) 

Jn[Z,  t)  —  DnCn  ^ 

(48) 

dxn(Z,t )  n  fd2Xn(Z,t)  |  1  3x„(Z,  t)\ 

3 1  ~Dn{  dz2  +z  dz  J 

(49) 

3 X„(z,t)  1  d[zjn(z,  t)] 

dt  Cn  Z  dz  ’ 

(50) 

wherejn(z,  t)  is  the  flux  of  hydrogen  diffusing  from  the  surface  to 
the  interior  of  the  electrode  at  spatial  position  z  and  time  t,  and 
Jn(zo ,  t)  =  I  ft  /(Fan).  If  each  cell  is  small  enough,  the  concentration 
x(z2,  t)inthei-  thcell(0  <  i  <  N),  (see  Figure(2)),  can  be  considered 
constant  with  input  and  output  hydrogen  flow  given  by  J(z2-,  t)  and 
J(zI+1 ,  t )  respectively.  Using  this  spatial  discretization  together  with 
the  forward  and  backward  approximation  of  the  derivative,  Eqs. 
(48)-(50)  can  be  written  as  follows: 

Uzh  0  =  -  ^  (x(Zi,  t)  -  x(Zj_i ,  t))  (51 ) 

^T* 1 11  =  “  e^Ai  ( WnUi+1 , 0  -  Zjnfr,  0)  ,  (52) 

Using  these  relationships  and  considering  a  boundary  condition  for 
the  flux  J(zN,  t)  =  0  at  the  active  material/metallic  substrate  inter¬ 
face  (z  =  Re),  the  diffusional  process  can  be  described  by  the  set  of 
Eqs.  (14)-(16). 


In  a  similar  way,  but  considering  the  spherical  geometry  of  the 
metal  hydride  particles,  the  equations  governing  the  H  concentra¬ 
tion  profile  in  the  MH  electrode  are  given  by: 

Jm(z,  t)  =  -Dm  Cm  dXm^z  f)  (53) 

dxm(z,t )  1  d[z2Jm(Z,t)} 

dt  cmz 2  dZ  ’  (  > 

where Jm (z,  t)Dm  is  the  diffusion  coefficient  and Jm(z0,  t)  =  If2/(Fam ) 
being  am  the  effective  transfer  area.  A  similar  discrete-spatial 
state  model  is  derived  for  the  metal  electrode  by  just  replacing 
A z  =  with  Rp  the  radium  of  the  sphere,  and  dt  =  z?+1  /z?  = 
(2i  +  3)2/(2i  +  l)2  for  1  <i<N. 
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